library(tidyverse) # tidy style coding
library(brms) # Bayesian models
library(bayesplot) # pretty bayes visuals
library(tidybayes) # Bayesian aesthetics
library(loo) # to use information criteria in brms models
library(MetBrewer) # colours
library(rcartocolor) # more colours
library(pander) # tables
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(ggdist) # for ribbon plot
library(png) # to load images
library(grid) # to plot images
\(~\)
Table S1. Recipe for food medium used in our experiment. The provided quantities make ~ 1 litre of food.
tibble("Ingredients" = c("Soy flour", "Cornmeal", "Yeast", "Dextrose", "Agar", "Water", "Tegosept", "Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL)"),
"Quantity" = c("20 g", "73 g", "35 g", "75 g", "6 g", "1000 mL", "17 mL", "14 mL")) %>%
pander(split.cell = 40, split.table = Inf)
| Ingredients | Quantity |
|---|---|
| Soy flour | 20 g |
| Cornmeal | 73 g |
| Yeast | 35 g |
| Dextrose | 75 g |
| Agar | 6 g |
| Water | 1000 mL |
| Tegosept | 17 mL |
| Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL) | 14 mL |
\(~\)
img <- readPNG("Figure_S1.png")
grid.raster(img)
Figure S1. Crossing scheme used to integrate the GFP constructs and apXA marked translocated second and third chromosome balancers into the LHM genetic background. We replicated the crosses 12 times to supply the flies used in generation zero of experimental evolution; 6 times using the Ubi GFP construct and 6 times with the 3xP construct. G = generation.
\(~\)
\(~\)
\(~\)
We conducted single-pair full-sibling crosses to erode heterozygosity across the genome, thereby exposing mutation load. Our proxies for mutation load - extinction rate and productivity decline - may also have been affected by differences in the intensity of interlocus sexual conflict between subpopulations carrying autosomes with different selection response histories. We specifically designed our experiment to minimise this effect, by enforcing monogamy between all breeding individuals across the extinction assay. However, we cannot discount males carrying male-limited autosomes and/ or control autosomes being more harmful to females than males carrying autosomes with a female-limited selection response history. This might occur because male harm is a by-product of male-male competition, which cannot respond to selection in the female-limited treatments. Furthermore, traits genetically correlated with harm such as activity levels or metabolic rate may be selected in the opposite direction among females, as these traits are highly sexually dimorphic. This may drive the evolution of reduced harm during female-limited evolution.
We observe 135 female (6.2%) and 21 (1%) male mortalities across the 2184 Drosophila vials used during the 20 generations of the experiment (we do not include the 9 final generations observed in Block 1 in this analysis). Male mortality did not occur often enough to provide the power required for formal analysis. We instead focus on modelling female mortality and test whether this is affected by selection response history.
We test this hypothesis to assess whether it is required to censor subpopulations where the breeding female died in the generation of extinction.
\(~\)
mortality_data <-
read_csv("Data/Mortality_data.csv") %>%
rowid_to_column("subpopulation") %>%
pivot_longer(cols = 8:27, names_to = "Generation", values_to = "Female_mortality") %>%
mutate(
Generation = as.integer(str_remove(Generation, "Gen_")),
Week = case_when(
Block == 1 ~ Generation,
Block == 2 ~ as.integer(Generation + 9)),
Week = as.factor(Week),
across(1:7, as.factor),
Male_mortality = Female_mortality,
Female_mortality = if_else(Female_mortality == "FEMALE" | Female_mortality == "BOTH", 1, 0),
Male_mortality = if_else(Male_mortality == "MALE" | Male_mortality == "BOTH", 1, 0)) %>%
rename(Evolution_treatment = Treatment)
# Create a function to build HTML searchable tables
my_data_table <- function(df){
datatable(
df, rownames=FALSE,
autoHideNavigation = TRUE,
extensions = c("Scroller", "Buttons"),
options = list(
dom = 'Bfrtip',
deferRender=TRUE,
scrollX=TRUE, scrollY=400,
scrollCollapse=TRUE,
buttons =
list('pageLength', 'colvis', 'csv', list(
extend = 'pdf',
pageSize = 'A4',
orientation = 'landscape',
filename = 'mortality_data')),
pageLength = 8500
)
)
}
my_data_table(mortality_data %>%
select(Mother_strain, Father_strain, Cross, subpopulation, Block, Week, Generation, Evolution_treatment, Female_mortality, Male_mortality))
\(~\)
Column explanations
Mother strain: we measured extinction and productivity
of families derived from 36 different crosses. Crosses were only
conducted between strains that experienced the same
Evolution_treatment. This column shows the strain that the
founding female of the subpopulation was derived from.
Father strain: we measured extinction and productivity
of families derived from 36 different crosses. Crosses were only
conducted between strains that experienced the same evolution treatment.
This column shows the strain that the founding male of the subpopulation
was derived from.
Cross: an identifying ID for each of the 36 combinations
of Mother strain and Father strain.
subpopulation: lines of flies founded by single
inseminated females in generation zero. We generated the next generation
of each subpopulation by crossing a single full-sibling dyad.
Block: the experiment was run in 2 distinct blocks,
using flies separated by 9 generations.
Week: the week that the generation of each block was
conducted in. Laboratory conditions may have differed slightly between
weeks. For example, a new batch of food was cooked for each week. Note
that blocks overlap between weeks 9-20.
Generation: the generation of the extinction assay,
ranging from 1 to 20.
Evolution_treatment: the strains had been exposed to one
of three evolutionary conditions for 20 generations: a female-limited
response to selection, a male-limited response and a control condition
where an evolutionary response occurred in both sexes.
Female_mortality: a 0 indicates the female survived
until she was removed from the 3-day breeding vial. A 1 indicates she
died during this period. Females were 1-4 days old when they entered the
breeding vial.
Male_mortality: see Female_mortality,
except this column documents male mortality in the breeding vials.
\(~\)
\(~\)
We fit a binomial model with Female_mortality as the
response, Evolution_treatment as a fixed effect and
Week, subpopulation and Cross as
random effects.
Priors
We fit moderately informative priors, with the aim to regularise our posterior estimates and improve model fitting by ruling out non-sensical values.
Note that each of these priors is expressed on the logit scale, due to the binomial logit link.
\(\alpha\) (the intercepts) ~ Normal(\(\mu\) = -2, \(\sigma\) = 2)
\(\beta\) (the fixed effects) ~ Normal(\(\mu\) = 0, \(\sigma\) = 2)
\(\sigma\) ~ Exponential(rate = 1)
female_mortality_model <-
brm(Female_mortality ~ 1 + Evolution_treatment + (1|Week) + (1|Cross) + (1|subpopulation),
data = mortality_data %>% filter(Female_mortality != "NA"),
family = bernoulli,
prior = c(prior(normal(-2, 2), class = Intercept),
prior(normal(0, 2), class = b),
prior(exponential(1), class = sd)),
iter = 8000, warmup = 4000, chains = 4, cores = 4, seed = 1,
control = list(adapt_delta = 0.98),
file = "Fits/female_mortality_model")
mortality_predictions <-
female_mortality_model %>%
as_draws_df() %>%
mutate(Control = inv_logit_scaled(b_Intercept)*100,
`Female-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentFemale)*100,
`Male-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentMale)*100) %>%
select(Control, `Female-limited`, `Male-limited`)
# put in easy to plot format
mortality_predictions_long <-
mortality_predictions %>%
pivot_longer(cols = everything(), names_to = "Evolution_treatment", values_to = "Mortality")
# calculate the differences between each treatment
mortality_diff <-
mortality_predictions %>%
mutate(`Control - Male-limited` = Control - `Male-limited`,
`Male-limited - Female-limited` = `Male-limited` - `Female-limited`,
`Control - Female-limited` = Control - `Female-limited`) %>%
pivot_longer(cols = 4:6, names_to = "Difference_contrast", values_to = "Difference") %>%
select(contains("Diff"))
mortality_plot <-
mortality_predictions_long %>%
ggplot(aes(x = Mortality, y = Evolution_treatment)) +
stat_halfeye(aes(fill = Evolution_treatment), .width = c(0.66, 0.95), alpha = 1,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
labs(x= "Female mortality (%)", y = "Inheritance treatment") +
theme_bw() +
theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
panel.grid.minor.x = element_blank(),
legend.position = "none", #transparent legend panel
text = element_text(size=14))
mortality_diff_plot <-
mortality_diff %>%
ggplot(aes(x = Difference, y = fct_relevel(Difference_contrast, "Control - Female-limited", "Male-limited - Female-limited", "Control - Male-limited"))) +
stat_halfeye(aes(fill = Difference_contrast), .width = c(0.66, 0.95), alpha = 1,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = c(carto_pal(7, "Peach")[2], carto_pal(7, "Purp")[1], carto_pal(7, "TealGrn")[1])) +
geom_vline(xintercept = 0, linetype = 2, colour = "black", linewidth = 1) +
labs(x = "Diff. in mortality (% points)", y = "Treatment contrast") +
scale_x_continuous(breaks=seq(-2, 6, 2), limits = c(-3, 7)) +
theme_bw() +
theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
panel.grid.minor.x = element_blank(),
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
legend.position = "none", #transparent legend panel
text = element_text(size=14))
mortality_plot / mortality_diff_plot +
plot_annotation(tag_levels = 'a')
Figure S2. Female mortality is less frequent in subpopulations with female-limited selection response histories, suggesting that male harm may be less intense in these subpopulations. Panel a shows the posterior distribution of the mean percentage of female mortality events across the total number of vials that housed subpopulations throughout the extinction assay, split by selection response history. Panel b shows the posterior distribution of the difference between each treatment. The points show the estimated median, with associated 66 and 95% credible intervals.
\(~\)
To avoid this confounding our measure of mutation load, we can censor extinction events that co-occur with female mortality.
\(~\)
\(~\)
extinction_data <-
read_csv("Data/Extinction_data.csv")
# Block 1 runs for 29 generations, while Block 2 only runs for 20. To calculate the censoring variable, we need to split these by Block, mutate the data, then rebind them
Block_1 <-
extinction_data %>%
filter(Block == "1") %>%
# here we create a censoring column. If the family a) escaped or was killed by something unrelated to the experiment or b) survived the 20 generations of the experiment, then we code a value of 1. If the family went extinct, we code a value of 0. This allows us to right censor the data, thereby preserving the information it provides on extinction.
mutate(across(Gen_1:Gen_29, ~replace_na(.x, "Escape")),
Censored_alive = if_else(Gen_29 == "YES", 1, 0),
Censored_escape = if_else(Gen_29 == "Escape", 1, 0),
Censored = Censored_alive + Censored_escape,
across(Gen_1:Gen_29, ~if_else(.x == "YES", 1, 0)))
Block_2 <-
extinction_data %>%
filter(Block == "2") %>%
# here we create a censoring column. If the family a) escaped or was killed by something unrelated to the experiment or b) survived the 20 generations of the experiment, then we code a value of 1. If the family went extinct, we code a value of 0. This allows us to right censor the data, thereby preserving the information it provides on extinction.
mutate(across(Gen_1:Gen_20, ~replace_na(.x, "Escape")),
across(Gen_21:Gen_29, ~replace_na(.x, "Not measured")),
Censored_alive = if_else(Gen_20 == "YES", 1, 0),
Censored_escape = if_else(Gen_20 == "Escape", 1, 0),
Censored = Censored_alive + Censored_escape,
across(Gen_1:Gen_29, ~if_else(.x == "YES", 1, 0)))
# combine the Blocked data back into a single tibble
extinction_data_wrangled <-
rbind(Block_1, Block_2) %>%
mutate(across(1:7, as.factor),
Gens_to_extinct = Gen_1 + Gen_2 + Gen_3 + Gen_4 +
Gen_5 + Gen_6 + Gen_7 + Gen_8 + Gen_9 + Gen_10 +
Gen_11 + Gen_12 + Gen_13 + Gen_14 + Gen_15 + Gen_16 +
Gen_17 + Gen_18 + Gen_19 + Gen_20 + Gen_21 + Gen_22 +
Gen_23 + Gen_24 + Gen_25 + Gen_26 + Gen_27 + Gen_28 + Gen_29 + 1) %>%
rename(Evolution_treatment = Treatment, subpopulation = ID) %>%
select(Mother_strain, Father_strain, Cross, subpopulation, Block,
Evolution_treatment, Gens_to_extinct, Gen_1:Gen_29, Censored_alive,
Censored_escape, Censored)
# Find the extinctions that co-occur with female mortality
extinction_mortality <-
left_join(
extinction_data_wrangled %>%
pivot_longer(cols = 8:36, names_to = "Generation", values_to = "Extant") %>%
mutate(Generation = as.integer(str_remove(Generation, "Gen_"))),
mortality_data
) %>%
filter(Extant == 0 & Female_mortality == 1) %>%
mutate(Censored_mortality = 1) %>%
select(subpopulation, Censored_mortality)
extinction_data_wrangled <- left_join(extinction_data_wrangled, extinction_mortality) %>%
mutate(Censored_mortality = if_else(is.na(Censored_mortality), 0, 1),
Censored_2 = Censored + Censored_mortality)
# Create a function to build HTML searchable tables
my_data_table <- function(df){
datatable(
df, rownames=FALSE,
autoHideNavigation = TRUE,
extensions = c("Scroller", "Buttons"),
options = list(
dom = 'Bfrtip',
deferRender=TRUE,
scrollX=TRUE, scrollY=400,
scrollCollapse=TRUE,
buttons =
list('pageLength', 'colvis', 'csv', list(
extend = 'pdf',
pageSize = 'A4',
orientation = 'landscape',
filename = 'extinction_data')),
pageLength = 500
)
)
}
my_data_table(extinction_data_wrangled)
Column explanations
Mother strain: we measured extinction and productivity
of families derived from 36 different crosses. Crosses were only
conducted between strains that experienced the same
Evolution_treatment. This column shows the strain that the
founding female of the subpopulation was derived from.
Father strain: we measured extinction and productivity
of families derived from 36 different crosses. Crosses were only
conducted between strains that experienced the same evolution treatment.
This column shows the strain that the founding male of the subpopulation
was derived from.
Cross: an identifying ID for each of the 36 combinations
of Mother strain and Father strain.
subpopulation: lines of flies founded by single
inseminated females in generation zero. We generated the next generation
of each subpopulation by crossing a single full-sibling dyad.
Block: the experiment was run in 2 distinct blocks,
using flies separated by 9 generations. Block 1 ran for 29 generations,
whereas Block 2 ran for 20 generations.
Evolution_treatment: the strains had been exposed to one
of three evolutionary conditions for 20-29 generations: a female-limited
response to selection, a male-limited response and a control condition
where an evolutionary response occurred in both sexes.
Gens_to_extinction: the number of generations the
subpopulation survived until it went extinct.
Gen_1:29: a 1 indicates the subpopulation was extant for
the generation in question, while a 0 indicates extinction. Note that we
continued propagating surviving subpopulations from Block 1 until
generation 29, which coincided with generation 20 for subpopulations
belonging to the second Block.
Censored: if the subpopulation 1) ended because flies
escaped or were killed by something unrelated to the experiment, 2) went
extinct because of breeding female mortality or 3) was extant after the
final generation of the experiment, then we coded a value of 1. If the
subpopulation went extinct, we coded a value of 0. This allows us to
right censor the data as per the brms syntax, thereby
preserving the information these censored cases provide on
extinction.
\(~\)
\(~\)
We fit a weibull survival model to estimate , the extinction rate of the subpopulations used in our experiment. The weibull model is an extension of an exponential decay model, and includes an additional shape parameter, which allows to vary across generations of inbreeding.
We right censor subpopulations that 1) ended because flies escaped or were killed by something unrelated to the experiment, 2) went extinct because of breeding female mortality or 3) were extant after the final generation of the experiment
Response variable
Gens_to_extinct is our response - a continuous variable
that runs from 1 to 30 generations. Note that we measured extinction up
until generation 29 in Block 1 and generation 20 in Block 2. A value of
30 indicates that the subpopulation was extant at the end of the
experiment in Block 1, while a value of 21 indicates this in Block 2. We
also include right censoring, to account for subpopulations extant at
the end of the experiment, or that were removed from the experiment by a
handling error. Censoring enables the inclusion of this ‘incomplete’
data in the model.
Fixed effects
Evolution_treatment: we include this to test for a
causal effect of sex-limited experimental evolution on mutation
load.
Block: extinction rate might differ between the blocks
we split our experiment up into e.g. because of microvariation in the
lab brought about because of a change of season, slight inconsistencies
in Drosophila food consistency etc.
Varying/Random effects
Cross: we initiated each subpopulation by crossing two
strains belonging to the same evolution treatment. There were 12 levels
of cross nested within evolution treatment, for a total of 36 crosses.
We kept subpopulations from the same Cross and
Block in the same column of the Drosophila vial
trays, so incuding this variable also controls for micro-environmental
variation caused by within-tray position. We used this random effect in
our first model.
Mother_strain and Father_strain: Used as an
alternative to Cross in our second model. In this case, we
fit a multi-membership model, where each subpopulation simultaneously
belongs to two levels of the Cross random effect. That is,
each subpopulation was founded my a mother from \(strain_i\) and a father from \(strain_j\) and they therefore always belong
to two levels of strain.
Multi-membership vs single-membership
While the multi-membership approach may initially seem the more powerful of the two random effect structures, we favour the single-membership approach for several reasons. First, our primary aim is to estimate the overall causal effect of evolution treatment on generations to extinction, rather than that of each strain (nested within evolution treatment). Given that we balanced our crossing design so that each strain is equally represented in the experiment, estimation for each strain is not integral. Second, Cross captures nuisance environmental variation that Strain does not. We distributed subpopulations throughout three Drosophila vial trays by placing subpopulations from the same cross in the same column. Across the columns we split subpopulations up by evolution treatment, such that column one contained subpopulations from the female-limited treatment, column two contained subpopulations from the male-limited treatment and column three contained subpopulations from the control treatment. We repeated this pattern until all 36 crosses filled a column. As a strain is used in multiple crosses, it is widely distirbuted throughout the trays and therefore does not capture a location effect like cross does. We are of the opinion that controlling for this environmental variation is of greater importance than strain level differences in mutation load, which both cross and evolution treatment also contain information for.
Finally, we fit both models and use leave one out (LOO) cross validation to formally assess which model has better expected predictive accuracy, both in and out of sample.
Priors
We fit moderately informative priors, with the aim to regularise our posterior estimates and improve model fitting by ruling out non-sensical values.
Note that each of these priors is expressed on the log scale, due to the weibull log link.
\(\alpha\) (the intercepts) ~ Gamma(shape = 5, rate = 3)
This intercept prior predicts an appropriate starting point for per generation extinction rate on the log scale.
\(\beta\) (the fixed effects) ~ Normal(\(\mu\) = 0, \(\sigma\) = 1)
\(\sigma\) ~ Exponential(rate = 1)
\(k\) (The weibull shape parameter) ~ Normal(\(\mu\) = 0, \(\sigma\) = 0.5)
extinction_model_sm <-
brm(data = extinction_data_wrangled,
family = weibull,
bf(Gens_to_extinct | cens(Censored_2) ~ 1 + Evolution_treatment + Block + (1|Cross),
shape ~ 1 + Evolution_treatment + Block),
# this intercept prior predicts an appropriate starting point for per gen extinction rate on the log scale
prior = c(prior(gamma(5, 3), class = Intercept, lb = 0),
prior(normal(0, 1), class = b),
prior(exponential(1), class = sd),
prior(normal(0, 0.5), class = b, dpar = shape)),
iter = 8000, warmup = 4000, chains = 4, cores = 4,
seed = 1, control = list(adapt_delta = .9, max_treedepth = 10),
file = "Fits/extinction_model")
extinction_model_sm <- add_criterion(extinction_model_sm, criterion = "loo")
# the model does not converge when we try and fit a treatment * block interaction
extinction_model_mm <-
brm(data = extinction_data_wrangled,
family = weibull,
bf(Gens_to_extinct | cens(Censored_2) ~ 1 + Evolution_treatment + Block + (1|mm(Mother_strain, Father_strain)),
shape ~ 1 + Evolution_treatment + Block),
# this intercept prior predicts an appropriate starting point for per gen extinction rate on the log scale
prior = c(prior(gamma(5, 3), class = Intercept, lb = 0),
prior(normal(0, 1), class = b),
prior(exponential(1), class = sd),
prior(normal(0, 0.5), class = b, dpar = shape)),
iter = 8000, warmup = 4000, chains = 4, cores = 4,
seed = 2, control = list(adapt_delta = 0.98, max_treedepth = 15),
file = "Fits/extinction_model_mm")
extinction_model_mm <- add_criterion(extinction_model_mm, criterion = "loo")
\(~\)
\(~\)
Compare models using LOO, a bayesian information
criteria
loo_compare(extinction_model_sm, extinction_model_mm) %>%
kable(digits = 3) %>%
kable_styling()
| elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | |
|---|---|---|---|---|---|---|---|---|
| extinction_model_sm | 0.000 | 0.000 | -1042.910 | 20.209 | 23.694 | 2.970 | 2085.820 | 40.418 |
| extinction_model_mm | -0.668 | 2.593 | -1043.578 | 20.238 | 15.072 | 2.273 | 2087.155 | 40.475 |
The single-membership model performs slightly better. We present results from this model.
Lets see how the model recapitulates the data
pp_check(extinction_model_sm, type = "hist", ndraws = 11, binwidth = 1) +
theme_minimal() +
theme(panel.background = element_blank())
\(~\)
\(~\)
We can initially ignore the shape parameter and find the
mean rate of extinction, averaged across all generations, using the
equation
\(\lambda= \frac{1}{e^\mu}\)
First find distributions of \(\mu\) for each of the three treatments.
inverse_rate_weibull <-
extinction_model_sm %>%
as_draws_df() %>%
transmute(`Control B1` = b_Intercept,
`Female_limited B1` = b_Intercept + b_Evolution_treatmentFemale,
`Male_limited B1` = b_Intercept + b_Evolution_treatmentMale,
`Control B2` = b_Intercept + b_Block2,
`Female_limited B2` = b_Intercept + b_Evolution_treatmentFemale + b_Block2,
`Male_limited B2` = b_Intercept + b_Evolution_treatmentMale + b_Block2)
Now lets plug these distributions into the equation to find mean estimates of \(\lambda\).
inverse_rate_weibull %>%
mutate(across(1:6, ~1 / exp(.x))) %>%
mutate(across(everything(), ~mean(.x))) %>%
distinct()
## # A tibble: 1 × 6
## `Control B1` `Female_limited B1` `Male_limited B1` Control B…¹ Femal…² Male_…³
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.142 0.130 0.111 0.120 0.110 0.0941
## # … with abbreviated variable names ¹`Control B2`, ²`Female_limited B2`,
## # ³`Male_limited B2`
This is the rate of extinction per generation, averaged across all generations. However, the weibull family allows a dynamic rate rather the the constant rate that the exponential constrains it to.
Following the vignette("brms_families"), the weibull
survival function can be found using the following transformation:
\(\lambda= \frac{e^\mu}{\gamma(1 + \frac{1}{k})}\)
where \(k\) is the shape parameter
estimated in brms
Then we can find survival at generation \(t\) using
\(S = e^{-(\frac{t}{\lambda})^k}\)
We plug in a vector with 16000 draws from the posterior distribution for 1) \(\mu\) of each inheritance treatment and 2) \(k\) the shape parameter to find distributions of the proportion of surviving subpopulations at each generation.
# Find lambda for the three treatments.
lambda_weibull <-
extinction_model_sm %>%
as_draws_df() %>%
transmute(`Control B1` = b_Intercept,
`Female_limited B1` = b_Intercept + b_Evolution_treatmentFemale,
`Male_limited B1` = b_Intercept + b_Evolution_treatmentMale,
`Control B2` = b_Intercept + b_Block2,
`Female_limited B2` = b_Intercept + b_Evolution_treatmentFemale + b_Block2,
`Male_limited B2` = b_Intercept + b_Evolution_treatmentMale + b_Block2,
shape_control_B1 = exp(b_shape_Intercept),
shape_female_B1 = exp(b_shape_Intercept + b_shape_Evolution_treatmentFemale),
shape_male_B1 = exp(b_shape_Intercept + b_shape_Evolution_treatmentMale),
shape_control_B2 = exp(b_shape_Intercept + b_shape_Block2),
shape_female_B2 = exp(b_shape_Intercept + b_shape_Evolution_treatmentFemale + b_shape_Block2),
shape_male_B2 = exp(b_shape_Intercept + b_shape_Evolution_treatmentMale + b_shape_Block2)) %>%
mutate(Control_lambda_B1 = exp(`Control B1`) / gamma(1 + 1/shape_control_B1),
Female_lambda_B1 = exp(`Female_limited B1`) / gamma(1 + 1/shape_female_B1),
Male_lambda_B1 = exp(`Male_limited B1`) / gamma(1 + 1/shape_male_B1),
Control_lambda_B2 = exp(`Control B2`) / gamma(1 + 1/shape_control_B2),
Female_lambda_B2 = exp(`Female_limited B2`) / gamma(1 + 1/shape_female_B2),
Male_lambda_B2 = exp(`Male_limited B2`) / gamma(1 + 1/shape_male_B2)) %>%
select(contains(c("lambda", "shape")))
# Now find survival proportion at generations 1-29
Weibull_survival_curve <-
lambda_weibull %>%
mutate(#Block 1
Control_0_B1 = exp(-(0/Control_lambda_B1)^shape_control_B1),
Control_1_B1 = exp(-(1/Control_lambda_B1)^shape_control_B1),
Control_2_B1 = exp(-(2/Control_lambda_B1)^shape_control_B1),
Control_3_B1 = exp(-(3/Control_lambda_B1)^shape_control_B1),
Control_4_B1 = exp(-(4/Control_lambda_B1)^shape_control_B1),
Control_5_B1 = exp(-(5/Control_lambda_B1)^shape_control_B1),
Control_6_B1 = exp(-(6/Control_lambda_B1)^shape_control_B1),
Control_7_B1 = exp(-(7/Control_lambda_B1)^shape_control_B1),
Control_8_B1 = exp(-(8/Control_lambda_B1)^shape_control_B1),
Control_9_B1 = exp(-(9/Control_lambda_B1)^shape_control_B1),
Control_10_B1 = exp(-(10/Control_lambda_B1)^shape_control_B1),
Control_11_B1 = exp(-(11/Control_lambda_B1)^shape_control_B1),
Control_12_B1 = exp(-(12/Control_lambda_B1)^shape_control_B1),
Control_13_B1 = exp(-(13/Control_lambda_B1)^shape_control_B1),
Control_14_B1 = exp(-(14/Control_lambda_B1)^shape_control_B1),
Control_15_B1 = exp(-(15/Control_lambda_B1)^shape_control_B1),
Control_16_B1 = exp(-(16/Control_lambda_B1)^shape_control_B1),
Control_17_B1 = exp(-(17/Control_lambda_B1)^shape_control_B1),
Control_18_B1 = exp(-(18/Control_lambda_B1)^shape_control_B1),
Control_19_B1 = exp(-(19/Control_lambda_B1)^shape_control_B1),
Control_20_B1 = exp(-(20/Control_lambda_B1)^shape_control_B1),
Control_21_B1 = exp(-(21/Control_lambda_B1)^shape_control_B1),
Control_22_B1 = exp(-(22/Control_lambda_B1)^shape_control_B1),
Control_23_B1 = exp(-(23/Control_lambda_B1)^shape_control_B1),
Control_24_B1 = exp(-(24/Control_lambda_B1)^shape_control_B1),
Control_25_B1 = exp(-(25/Control_lambda_B1)^shape_control_B1),
Control_26_B1 = exp(-(26/Control_lambda_B1)^shape_control_B1),
Control_27_B1 = exp(-(27/Control_lambda_B1)^shape_control_B1),
Control_28_B1 = exp(-(28/Control_lambda_B1)^shape_control_B1),
Control_29_B1 = exp(-(29/Control_lambda_B1)^shape_control_B1),
Female_0_B1 = exp(-(0/Female_lambda_B1)^shape_female_B1),
Female_1_B1 = exp(-(1/Female_lambda_B1)^shape_female_B1),
Female_2_B1 = exp(-(2/Female_lambda_B1)^shape_female_B1),
Female_3_B1 = exp(-(3/Female_lambda_B1)^shape_female_B1),
Female_4_B1 = exp(-(4/Female_lambda_B1)^shape_female_B1),
Female_5_B1 = exp(-(5/Female_lambda_B1)^shape_female_B1),
Female_6_B1 = exp(-(6/Female_lambda_B1)^shape_female_B1),
Female_7_B1 = exp(-(7/Female_lambda_B1)^shape_female_B1),
Female_8_B1 = exp(-(8/Female_lambda_B1)^shape_female_B1),
Female_9_B1 = exp(-(9/Female_lambda_B1)^shape_female_B1),
Female_10_B1 = exp(-(10/Female_lambda_B1)^shape_female_B1),
Female_11_B1 = exp(-(11/Female_lambda_B1)^shape_female_B1),
Female_12_B1 = exp(-(12/Female_lambda_B1)^shape_female_B1),
Female_13_B1 = exp(-(13/Female_lambda_B1)^shape_female_B1),
Female_14_B1 = exp(-(14/Female_lambda_B1)^shape_female_B1),
Female_15_B1 = exp(-(15/Female_lambda_B1)^shape_female_B1),
Female_16_B1 = exp(-(16/Female_lambda_B1)^shape_female_B1),
Female_17_B1 = exp(-(17/Female_lambda_B1)^shape_female_B1),
Female_18_B1 = exp(-(18/Female_lambda_B1)^shape_female_B1),
Female_19_B1 = exp(-(19/Female_lambda_B1)^shape_female_B1),
Female_20_B1 = exp(-(20/Female_lambda_B1)^shape_female_B1),
Female_21_B1 = exp(-(21/Female_lambda_B1)^shape_female_B1),
Female_22_B1 = exp(-(22/Female_lambda_B1)^shape_female_B1),
Female_23_B1 = exp(-(23/Female_lambda_B1)^shape_female_B1),
Female_24_B1 = exp(-(24/Female_lambda_B1)^shape_female_B1),
Female_25_B1 = exp(-(25/Female_lambda_B1)^shape_female_B1),
Female_26_B1 = exp(-(26/Female_lambda_B1)^shape_female_B1),
Female_27_B1 = exp(-(27/Female_lambda_B1)^shape_female_B1),
Female_28_B1 = exp(-(28/Female_lambda_B1)^shape_female_B1),
Female_29_B1 = exp(-(29/Female_lambda_B1)^shape_female_B1),
Male_0_B1 = exp(-(0/Male_lambda_B1)^shape_male_B1),
Male_1_B1 = exp(-(1/Male_lambda_B1)^shape_male_B1),
Male_2_B1 = exp(-(2/Male_lambda_B1)^shape_male_B1),
Male_3_B1 = exp(-(3/Male_lambda_B1)^shape_male_B1),
Male_4_B1 = exp(-(4/Male_lambda_B1)^shape_male_B1),
Male_5_B1 = exp(-(5/Male_lambda_B1)^shape_male_B1),
Male_6_B1 = exp(-(6/Male_lambda_B1)^shape_male_B1),
Male_7_B1 = exp(-(7/Male_lambda_B1)^shape_male_B1),
Male_8_B1 = exp(-(8/Male_lambda_B1)^shape_male_B1),
Male_9_B1 = exp(-(9/Male_lambda_B1)^shape_male_B1),
Male_10_B1 = exp(-(10/Male_lambda_B1)^shape_male_B1),
Male_11_B1 = exp(-(11/Male_lambda_B1)^shape_male_B1),
Male_12_B1 = exp(-(12/Male_lambda_B1)^shape_male_B1),
Male_13_B1 = exp(-(13/Male_lambda_B1)^shape_male_B1),
Male_14_B1 = exp(-(14/Male_lambda_B1)^shape_male_B1),
Male_15_B1 = exp(-(15/Male_lambda_B1)^shape_male_B1),
Male_16_B1 = exp(-(16/Male_lambda_B1)^shape_male_B1),
Male_17_B1 = exp(-(17/Male_lambda_B1)^shape_male_B1),
Male_18_B1 = exp(-(18/Male_lambda_B1)^shape_male_B1),
Male_19_B1 = exp(-(19/Male_lambda_B1)^shape_male_B1),
Male_20_B1 = exp(-(20/Male_lambda_B1)^shape_male_B1),
Male_21_B1 = exp(-(21/Male_lambda_B1)^shape_male_B1),
Male_22_B1 = exp(-(22/Male_lambda_B1)^shape_male_B1),
Male_23_B1 = exp(-(23/Male_lambda_B1)^shape_male_B1),
Male_24_B1 = exp(-(24/Male_lambda_B1)^shape_male_B1),
Male_25_B1 = exp(-(25/Male_lambda_B1)^shape_male_B1),
Male_26_B1 = exp(-(26/Male_lambda_B1)^shape_male_B1),
Male_27_B1 = exp(-(27/Male_lambda_B1)^shape_male_B1),
Male_28_B1 = exp(-(28/Male_lambda_B1)^shape_male_B1),
Male_29_B1 = exp(-(29/Male_lambda_B1)^shape_male_B1),
# Block 2
Control_0_B2 = exp(-(0/Control_lambda_B2)^shape_control_B2),
Control_1_B2 = exp(-(1/Control_lambda_B2)^shape_control_B2),
Control_2_B2 = exp(-(2/Control_lambda_B2)^shape_control_B2),
Control_3_B2 = exp(-(3/Control_lambda_B2)^shape_control_B2),
Control_4_B2 = exp(-(4/Control_lambda_B2)^shape_control_B2),
Control_5_B2 = exp(-(5/Control_lambda_B2)^shape_control_B2),
Control_6_B2 = exp(-(6/Control_lambda_B2)^shape_control_B2),
Control_7_B2 = exp(-(7/Control_lambda_B2)^shape_control_B2),
Control_8_B2 = exp(-(8/Control_lambda_B2)^shape_control_B2),
Control_9_B2 = exp(-(9/Control_lambda_B2)^shape_control_B2),
Control_10_B2 = exp(-(10/Control_lambda_B2)^shape_control_B2),
Control_11_B2 = exp(-(11/Control_lambda_B2)^shape_control_B2),
Control_12_B2 = exp(-(12/Control_lambda_B2)^shape_control_B2),
Control_13_B2 = exp(-(13/Control_lambda_B2)^shape_control_B2),
Control_14_B2 = exp(-(14/Control_lambda_B2)^shape_control_B2),
Control_15_B2 = exp(-(15/Control_lambda_B2)^shape_control_B2),
Control_16_B2 = exp(-(16/Control_lambda_B2)^shape_control_B2),
Control_17_B2 = exp(-(17/Control_lambda_B2)^shape_control_B2),
Control_18_B2 = exp(-(18/Control_lambda_B2)^shape_control_B2),
Control_19_B2 = exp(-(19/Control_lambda_B2)^shape_control_B2),
Control_20_B2 = exp(-(20/Control_lambda_B2)^shape_control_B2),
Control_21_B2 = exp(-(21/Control_lambda_B2)^shape_control_B2),
Control_22_B2 = exp(-(22/Control_lambda_B2)^shape_control_B2),
Control_23_B2 = exp(-(23/Control_lambda_B2)^shape_control_B2),
Control_24_B2 = exp(-(24/Control_lambda_B2)^shape_control_B2),
Control_25_B2 = exp(-(25/Control_lambda_B2)^shape_control_B2),
Control_26_B2 = exp(-(26/Control_lambda_B2)^shape_control_B2),
Control_27_B2 = exp(-(27/Control_lambda_B2)^shape_control_B2),
Control_28_B2 = exp(-(28/Control_lambda_B2)^shape_control_B2),
Control_29_B2 = exp(-(29/Control_lambda_B2)^shape_control_B2),
Female_0_B2 = exp(-(0/Female_lambda_B2)^shape_female_B2),
Female_1_B2 = exp(-(1/Female_lambda_B2)^shape_female_B2),
Female_2_B2 = exp(-(2/Female_lambda_B2)^shape_female_B2),
Female_3_B2 = exp(-(3/Female_lambda_B2)^shape_female_B2),
Female_4_B2 = exp(-(4/Female_lambda_B2)^shape_female_B2),
Female_5_B2 = exp(-(5/Female_lambda_B2)^shape_female_B2),
Female_6_B2 = exp(-(6/Female_lambda_B2)^shape_female_B2),
Female_7_B2 = exp(-(7/Female_lambda_B2)^shape_female_B2),
Female_8_B2 = exp(-(8/Female_lambda_B2)^shape_female_B2),
Female_9_B2 = exp(-(9/Female_lambda_B2)^shape_female_B2),
Female_10_B2 = exp(-(10/Female_lambda_B2)^shape_female_B2),
Female_11_B2 = exp(-(11/Female_lambda_B2)^shape_female_B2),
Female_12_B2 = exp(-(12/Female_lambda_B2)^shape_female_B2),
Female_13_B2 = exp(-(13/Female_lambda_B2)^shape_female_B2),
Female_14_B2 = exp(-(14/Female_lambda_B2)^shape_female_B2),
Female_15_B2 = exp(-(15/Female_lambda_B2)^shape_female_B2),
Female_16_B2 = exp(-(16/Female_lambda_B2)^shape_female_B2),
Female_17_B2 = exp(-(17/Female_lambda_B2)^shape_female_B2),
Female_18_B2 = exp(-(18/Female_lambda_B2)^shape_female_B2),
Female_19_B2 = exp(-(19/Female_lambda_B2)^shape_female_B2),
Female_20_B2 = exp(-(20/Female_lambda_B2)^shape_female_B2),
Female_21_B2 = exp(-(21/Female_lambda_B2)^shape_female_B2),
Female_22_B2 = exp(-(22/Female_lambda_B2)^shape_female_B2),
Female_23_B2 = exp(-(23/Female_lambda_B2)^shape_female_B2),
Female_24_B2 = exp(-(24/Female_lambda_B2)^shape_female_B2),
Female_25_B2 = exp(-(25/Female_lambda_B2)^shape_female_B2),
Female_26_B2 = exp(-(26/Female_lambda_B2)^shape_female_B2),
Female_27_B2 = exp(-(27/Female_lambda_B2)^shape_female_B2),
Female_28_B2 = exp(-(28/Female_lambda_B2)^shape_female_B2),
Female_29_B2 = exp(-(29/Female_lambda_B2)^shape_female_B2),
Male_0_B2 = exp(-(0/Male_lambda_B2)^shape_male_B2),
Male_1_B2 = exp(-(1/Male_lambda_B2)^shape_male_B2),
Male_2_B2 = exp(-(2/Male_lambda_B2)^shape_male_B2),
Male_3_B2 = exp(-(3/Male_lambda_B2)^shape_male_B2),
Male_4_B2 = exp(-(4/Male_lambda_B2)^shape_male_B2),
Male_5_B2 = exp(-(5/Male_lambda_B2)^shape_male_B2),
Male_6_B2 = exp(-(6/Male_lambda_B2)^shape_male_B2),
Male_7_B2 = exp(-(7/Male_lambda_B2)^shape_male_B2),
Male_8_B2 = exp(-(8/Male_lambda_B2)^shape_male_B2),
Male_9_B2 = exp(-(9/Male_lambda_B2)^shape_male_B2),
Male_10_B2 = exp(-(10/Male_lambda_B2)^shape_male_B2),
Male_11_B2 = exp(-(11/Male_lambda_B2)^shape_male_B2),
Male_12_B2 = exp(-(12/Male_lambda_B2)^shape_male_B2),
Male_13_B2 = exp(-(13/Male_lambda_B2)^shape_male_B2),
Male_14_B2 = exp(-(14/Male_lambda_B2)^shape_male_B2),
Male_15_B2 = exp(-(15/Male_lambda_B2)^shape_male_B2),
Male_16_B2 = exp(-(16/Male_lambda_B2)^shape_male_B2),
Male_17_B2 = exp(-(17/Male_lambda_B2)^shape_male_B2),
Male_18_B2 = exp(-(18/Male_lambda_B2)^shape_male_B2),
Male_19_B2 = exp(-(19/Male_lambda_B2)^shape_male_B2),
Male_20_B2 = exp(-(20/Male_lambda_B2)^shape_male_B2),
Male_21_B2 = exp(-(21/Male_lambda_B2)^shape_male_B2),
Male_22_B2 = exp(-(22/Male_lambda_B2)^shape_male_B2),
Male_23_B2 = exp(-(23/Male_lambda_B2)^shape_male_B2),
Male_24_B2 = exp(-(24/Male_lambda_B2)^shape_male_B2),
Male_25_B2 = exp(-(25/Male_lambda_B2)^shape_male_B2),
Male_26_B2 = exp(-(26/Male_lambda_B2)^shape_male_B2),
Male_27_B2 = exp(-(27/Male_lambda_B2)^shape_male_B2),
Male_28_B2 = exp(-(28/Male_lambda_B2)^shape_male_B2),
Male_29_B2 = exp(-(29/Male_lambda_B2)^shape_male_B2)) %>%
select(-c(contains(c("shape", "lambda"))))
Weibull_extinction_estimates_long <-
Weibull_survival_curve %>%
pivot_longer(cols = everything(), names_to = "Evolution treatment", values_to = "Prop_subpopulations_surviving") %>%
separate(`Evolution treatment`, into = c("Evolution treatment", "Generation", "Block"), sep = "_") %>%
mutate(`Evolution treatment` = case_when(
`Evolution treatment` == "Female" ~ "Female-limited",
`Evolution treatment` == "Male" ~ "Male-limited",
`Evolution treatment` == "Control" ~ "Control"
))
\(~\)
Make panel a of Figure 1
# values from Falconer and Mackay 1996 pg 90
Inbreeding_coefs <-
tibble(Generation = 0:20,
Probability_homozygous = c(0, 0.25, 0.375, 0.5, 0.594, 0.672, 0.734, 0.785, 0.826,
0.859, 0.886, 0.908, 0.926, 0.94, 0.951, 0.961, 0.968, 0.974,
0.979, 0.983, 0.986),
Probability_pop_fixation = c(0, 0, 0.063, 0.172, 0.293, 0.409, 0.512, 0.601, 0.675, 0.736, 0.785,
0.826, 0.859, 0.886, 0.908, 0.925, 0.94, 0.951, 0.96, 0.968, 0.975)) %>%
mutate(single_locus_heterozygosity = 1 - Probability_homozygous,
single_locus_not_fixed = 1 - Probability_pop_fixation)
# for the plot annotations
arrows <-
tibble(
x1 = 3.2,
x2 = 5.6,
y1 = 0.28,
y2 = 0.49)
# Block 1
weibull_surv_plot_B1 <-
Weibull_extinction_estimates_long %>%
filter(Block == "B1") %>%
group_by(`Evolution treatment`, Generation) %>%
tidybayes::median_qi(Prop_subpopulations_surviving, .width = 0.5) %>%
mutate(Generation = as.numeric(Generation)) %>%
arrange(Generation) %>%
# plot!
ggplot(aes(x = Generation)) +
geom_line(data = Inbreeding_coefs, aes(y = single_locus_not_fixed), linetype = 3,
linewidth = 2, alpha = 1) +
geom_ribbon(aes(ymin = .lower, ymax = .upper, fill = `Evolution treatment`),
alpha = 1/2) +
geom_line(aes(y = Prop_subpopulations_surviving, color = `Evolution treatment`), linetype =5, linewidth = 0.8) +
annotate("text", x = 3.2, y = 0.26, size = 4, label = "Single locus prob. of heterozygosity") +
geom_curve(
data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
arrow = arrow(length = unit(0.08, "inch")), size = 0.75,
color = "gray20", curvature = -0.3) +
scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
scale_color_manual(values = met.brewer("Hiroshige", 3), guide = "none") +
scale_x_continuous(breaks = c(0, 5, 10, 15, 20, 25, 30), expand = expansion(mult = c(0.01, 0.01)), limits = c(0, 20)) +
scale_y_continuous(breaks = c(0, 0.25, .5, 0.75, 1), limits = 0:1, expand = expansion(mult = c(0.01, 0.01))) +
labs(x = "Generations of inbreeding", y = "Proportion of surviving subpopulations", fill = "Inheritance treatment") +
theme_bw() +
theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
legend.background = element_rect(fill='transparent'), #transparent legend bg
legend.box.background = element_rect(fill='transparent', colour = 'transparent'), #transparent legend panel
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6),
text = element_text(size = 14),
legend.title = element_text(size = 12),
legend.text = element_text(size = 11))
We can also plot the mean generations until extinction for each evolution treatment. Some data wrangling is first needed.
new_data <- expand_grid(
Evolution_treatment = extinction_data_wrangled$Evolution_treatment,
Block = 1:2) %>%
distinct(Evolution_treatment, Block) %>%
mutate(key = paste("V", 1:n(), sep = ""))
weibull_estimates <-
fitted(extinction_model_sm, newdata = new_data,
re_formula = NA, summary = F) %>% #robust = T) %>%
as_tibble() %>%
rename(`Female-limited_1` = V1, `Female-limited_2` = V2, `Male-limited_1` = V3, `Male-limited_2` = V4, Control_1 = V5, Control_2 = V6)
mean_weibull_estimates <-
weibull_estimates %>%
pivot_longer(cols = 1:6, names_to = "Evolution treatment", values_to = "generations") %>%
separate(col = `Evolution treatment`, into = c("Evolution treatment", "Block"), sep = "_")
# calculate contrasts
weibull_contrasts_B1 <-
weibull_estimates %>%
mutate(`Female - Control` = `Female-limited_1` - Control_1,
`Male - Female` = `Male-limited_1` - `Female-limited_1`,
`Male - Control` = `Male-limited_1` - Control_1) %>%
select(-c(`Male-limited_1`, Control_1, `Female-limited_1`,
`Male-limited_2`, Control_2, `Female-limited_2`)) %>%
pivot_longer(cols = 1:3, names_to = "Contrast", values_to = "generation_diff")
weibull_contrasts_B2 <-
weibull_estimates %>%
mutate(`Female - Control` = `Female-limited_2` - Control_2,
`Male - Female` = `Male-limited_2` - `Female-limited_2`,
`Male - Control` = `Male-limited_2` - Control_2) %>%
select(-c(`Male-limited_1`, Control_1, `Female-limited_1`,
`Male-limited_2`, Control_2, `Female-limited_2`)) %>%
pivot_longer(cols = 1:3, names_to = "Contrast", values_to = "generation_diff")
Create panels b and c of Figure 1
# plot the means
# block 1
ext_p1_B1 <-
mean_weibull_estimates %>%
filter(Block == "1") %>%
ggplot(aes(x = generations, y = `Evolution treatment`)) +
stat_halfeye(aes(fill = `Evolution treatment`), .width = c(0.66, 0.95), alpha = 1,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals + # width indicates the uncertainty intervals: here we have 66% and 95% intervals+
scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
labs(x = "Generations till extinction", y = "Inheritance treatment") +
scale_x_continuous(breaks=seq(4, 12, 1), limits = c(4, 12)) +
theme_bw() +
theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
panel.grid.minor.x = element_blank(),
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
legend.position = "none", #transparent legend panel
text = element_text(size=14))
# plot the contrasts
ext_p2_B1 <-
weibull_contrasts_B1 %>%
ggplot(aes(x = generation_diff, y = fct_relevel(Contrast, "Female - Control", "Male - Female", "Male - Control"))) +
stat_halfeye(aes(fill = Contrast), .width = c(0.66, 0.95), alpha = 1,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = c(carto_pal(7, "Peach")[2], carto_pal(7, "Purp")[1], carto_pal(7, "TealGrn")[1])) +
geom_vline(xintercept = 0, linetype = 2, colour = "black", size = 1) +
labs(x = "Diff. in generations till extinction", y = "Treatment contrast") +
scale_x_continuous(breaks=seq(-4, 6, 2)) +
theme_bw() +
theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
panel.grid.minor.x = element_blank(),
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
legend.position = "none", #transparent legend panel
text = element_text(size=14))
We can also estimate the difference between evolutionary treatments in the proportion of subpopulations extant in each generation. Below we calculate each contrast and make panels d, e and f
# Block 1
Difference_gens_B1 <-
Weibull_survival_curve %>%
select(contains("B1")) %>%
mutate(G1_diff.c = Male_1_B1 - Control_1_B1,
G1_diff.f = Male_1_B1 - Female_1_B1,
G1_diff.fc = Female_1_B1 - Control_1_B1,
G2_diff.c = Male_2_B1 - Control_2_B1,
G2_diff.f = Male_2_B1 - Female_2_B1,
G2_diff.fc = Female_2_B1 - Control_2_B1,
G3_diff.c = Male_3_B1 - Control_3_B1,
G3_diff.f = Male_3_B1 - Female_3_B1,
G3_diff.fc = Female_3_B1 - Control_3_B1,
G4_diff.c = Male_4_B1 - Control_4_B1,
G4_diff.f = Male_4_B1 - Female_4_B1,
G4_diff.fc = Female_4_B1 - Control_4_B1,
G5_diff.c = Male_5_B1 - Control_5_B1,
G5_diff.f = Male_5_B1 - Female_5_B1,
G5_diff.fc = Female_5_B1 - Control_5_B1,
G6_diff.c = Male_6_B1 - Control_6_B1,
G6_diff.f = Male_6_B1 - Female_6_B1,
G6_diff.fc = Female_6_B1 - Control_6_B1,
G7_diff.c = Male_7_B1 - Control_7_B1,
G7_diff.f = Male_7_B1 - Female_7_B1,
G7_diff.fc = Female_7_B1 - Control_7_B1,
G8_diff.c = Male_8_B1 - Control_8_B1,
G8_diff.f = Male_8_B1 - Female_8_B1,
G8_diff.fc = Female_8_B1 - Control_8_B1,
G9_diff.c = Male_9_B1 - Control_9_B1,
G9_diff.f = Male_9_B1 - Female_9_B1,
G9_diff.fc = Female_9_B1 - Control_9_B1,
G10_diff.c = Male_10_B1 - Control_10_B1,
G10_diff.f = Male_10_B1 - Female_10_B1,
G10_diff.fc = Female_10_B1 - Control_10_B1,
G11_diff.c = Male_11_B1 - Control_11_B1,
G11_diff.f = Male_11_B1 - Female_11_B1,
G11_diff.fc = Female_11_B1 - Control_11_B1,
G12_diff.c = Male_12_B1 - Control_12_B1,
G12_diff.f = Male_12_B1 - Female_12_B1,
G12_diff.fc = Female_12_B1 - Control_12_B1,
G13_diff.c = Male_13_B1 - Control_13_B1,
G13_diff.f = Male_13_B1 - Female_13_B1,
G13_diff.fc = Female_13_B1 - Control_13_B1,
G14_diff.c = Male_14_B1 - Control_14_B1,
G14_diff.f = Male_14_B1 - Female_14_B1,
G14_diff.fc = Female_14_B1 - Control_14_B1,
G15_diff.c = Male_15_B1 - Control_15_B1,
G15_diff.f = Male_15_B1 - Female_15_B1,
G15_diff.fc = Female_15_B1 - Control_15_B1,
G16_diff.c = Male_16_B1 - Control_16_B1,
G16_diff.f = Male_16_B1 - Female_16_B1,
G16_diff.fc = Female_16_B1 - Control_16_B1,
G17_diff.c = Male_17_B1 - Control_17_B1,
G17_diff.f = Male_17_B1 - Female_17_B1,
G17_diff.fc = Female_17_B1 - Control_17_B1,
G18_diff.c = Male_18_B1 - Control_18_B1,
G18_diff.f = Male_18_B1 - Female_18_B1,
G18_diff.fc = Female_18_B1 - Control_18_B1,
G19_diff.c = Male_19_B1 - Control_19_B1,
G19_diff.f = Male_19_B1 - Female_19_B1,
G19_diff.fc = Female_19_B1 - Control_19_B1,
G20_diff.c = Male_20_B1 - Control_20_B1,
G20_diff.f = Male_20_B1 - Female_20_B1,
G20_diff.fc = Female_20_B1 - Control_20_B1) %>%
#G21_diff.c = Male_21_B1 - Control_21_B1,
#G21_diff.f = Male_21_B1 - Female_21_B1,
#G21_diff.fc = Female_21_B1 - Control_21_B1,
#G22_diff.c = Male_22_B1 - Control_22_B1,
#G22_diff.f = Male_22_B1 - Female_22_B1,
#G22_diff.fc = Female_22_B1 - Control_22_B1,
#G23_diff.c = Male_23_B1 - Control_23_B1,
#G23_diff.f = Male_23_B1 - Female_23_B1,
#G23_diff.fc = Female_23_B1 - Control_23_B1,
#G24_diff.c = Male_24_B1 - Control_24_B1,
#G24_diff.f = Male_24_B1 - Female_24_B1,
#G24_diff.fc = Female_24_B1 - Control_24_B1,
#G25_diff.c = Male_25_B1 - Control_25_B1,
#G25_diff.f = Male_25_B1 - Female_25_B1,
#G25_diff.fc = Female_25_B1 - Control_25_B1,
#G26_diff.c = Male_26_B1 - Control_26_B1,
#G26_diff.f = Male_26_B1 - Female_26_B1,
#G26_diff.fc = Female_26_B1 - Control_26_B1,
#G27_diff.c = Male_27_B1 - Control_27_B1,
#G27_diff.f = Male_27_B1 - Female_27_B1,
#G27_diff.fc = Female_27_B1 - Control_27_B1,
#G28_diff.c = Male_28_B1 - Control_28_B1,
#G28_diff.f = Male_28_B1 - Female_28_B1,
#G28_diff.fc = Female_28_B1 - Control_28_B1,
#G29_diff.c = Male_2_B1 - Control_29_B1,
#G29_diff.f = Male_29_B1 - Female_29_B1,
#G29_diff.fc = Female_29_B1 - Control_29_B1) %>%
select(starts_with("G")) %>%
pivot_longer(cols = everything(), names_to = "Difference contrast", values_to = "survival_diff") %>%
separate(`Difference contrast`, into = c("Generation", "Difference contrast"), sep = "_") %>%
mutate(Generation = as.numeric(str_remove(Generation, "G")))
# Make the three plots
Leaf_plot_mc_B1 <-
Difference_gens_B1 %>%
filter(`Difference contrast` == "diff.c") %>%
ggplot(aes(survival_diff, as.factor(Generation))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Purp") +
coord_cartesian(xlim = c(-0.1, 0.3)) +
scale_y_discrete(limits=rev) +
geom_vline(linetype = 2, xintercept = 0, size = 1) +
labs(x = "Male - Control\nProp. surv subpopulations",
y = "Generation of inbreeding") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
Leaf_plot_mf_B1 <-
Difference_gens_B1 %>%
filter(`Difference contrast` == "diff.f") %>%
ggplot(aes(survival_diff, as.factor(Generation))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "TealGrn") +
coord_cartesian(xlim = c(-0.1, 0.3)) +
scale_y_discrete(limits=rev) +
geom_vline(linetype = 2, xintercept = 0, size = 1) +
labs(x = "Male - Female\nProp. surv subpopulations",
y = NULL) +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
Leaf_plot_fc_B1 <-
Difference_gens_B1 %>%
filter(`Difference contrast` == "diff.fc") %>%
ggplot(aes(survival_diff, as.factor(Generation))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Peach") +
coord_cartesian(xlim = c(-0.1, 0.3)) +
scale_y_discrete(limits=rev) +
geom_vline(linetype = 2, xintercept = 0, size = 1) +
labs(x = "Female - Control\nProp. surv subpopulations",
y = NULL) +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
Combine the 6 panels
(weibull_surv_plot_B1 / (ext_p1_B1 + ext_p2_B1)) / (Leaf_plot_mc_B1 + Leaf_plot_mf_B1 + Leaf_plot_fc_B1) +
plot_annotation(tag_levels = 'a') +
plot_layout(heights = c(1.5, 1, 1))
Figure 1. Panel a shows extinction
trajectories under increasing stress from inbreeding. Bands are 50%
credible intervals to aid visual inspection. The black dotted line is
the probability that a randomly selected autosomal locus remains
heterozygous. b shows the posterior distribution for
the number of generations that sub-populations descended from each
inheritance treatment survived inbreeding exposure. The white point
shows the median with associated 66 and 95% credible intervals c shows
the distribution of the difference between each treatment in generations
survived. d-f show the distribution of the difference
between each treatment in generations survived, estimated separately for
each generation of inbreeding. The innermost colour band approximates
the median, while the middle and outermost bars show 66 and 95% credible
intervals.
\(~\)
\(~\)
We have several hypotheses that we test with the productivity data. First, productivity should be negatively correlated with the number of generations inbreeding occurs, due to inbreeding depression. As mentioned above, this effect should be strongest in the first ~10 generations of the experiment, during which full-sibling inbreeding should have the largest effects on genome wide heterozygosity. After this point, the inbreeding coefficient will be 0.89, while over the next 10 generations it should only increase by a further 0.1 to 0.99. Hence, we expect productivity to stabilise between the 10th and 20th generations, assuming the subpopulation is still extant. Second, subpopulations carrying autosomes that have responded to selection on males should exhibit a smaller drop in productivity compared with subpopulations carrying the female-adapted or control autosomes, because they have a history of stronger selection that should, in theory, have purged the genome of recessive deleterious alleles.
data <-
read_csv("Data/Productivity_data.csv") %>%
mutate(Week = case_when(
Block == 1 ~ Generation,
Block == 2 ~ Generation + 9))
Productivity_data <-
left_join(
data,
data %>%
distinct(Cross, Replicate) %>%
rowid_to_column("subpopulation")
) %>%
select(subpopulation, everything()) %>%
mutate(across(1:7, as.factor),
Week = as.factor(Week),
Collection_window_offspring = Female_offspring + Male_offspring,
Pre_window_offspring = Pre_window_female_offspring + Pre_window_male_offspring,
Total_female_offspring = Female_offspring + Pre_window_female_offspring,
Total_male_offspring = Male_offspring + Pre_window_male_offspring,
Total_offspring = Total_female_offspring + Total_male_offspring) %>%
# In one generation of the experiment (b1 = G24 & B2 = G15) sibling pairs were setup a day early and removed from their vials at the regular time, meaning that they had an extra day to produce offspring. To correct for this we multiply offspring counts by 0.75.
mutate(Collection_window_offspring = if_else(Count_conditions == "Extra day",
round(Collection_window_offspring * 0.75), Collection_window_offspring),
# note that because everything is moved a day early, pre-window offspring counts will still be inflated even after this correction
Pre_window_offspring = if_else(Count_conditions == "Extra day",
round(Pre_window_offspring * 0.75), Pre_window_offspring))
# Combine with mortality and extinction data to filter out vials that went extinct because of female mortality
Productivity_mortality <-
left_join(
Productivity_data %>%
filter(Generation < 21), #%>% only include data collected on the first 20 generations of inbreeding (block 2's endpoint)
mortality_data %>% select(subpopulation, Generation, Female_mortality)
)
# We also include a column that specifies if a subpopulation was extinct. This means we can easily remove 0 values from the data if required. We can use the `extinction_data_wrangled$Gens_to_extinct` column to help us here.
Productivity_data_clean <-
left_join(
Productivity_mortality,
extinction_data_wrangled %>%
select(subpopulation, Gens_to_extinct, Censored_mortality)
) %>%
mutate(Vial_setup = if_else(Generation > Gens_to_extinct, "No", "Yes")) %>%
filter(Collection_window_offspring != "NA")
# remove mortality events and their consequences in later generations
Productivity_data_clean_2 <-
Productivity_data_clean %>%
mutate(Collection_window_offspring = case_when(
Censored_mortality == 1 & Vial_setup == "No" ~ NA_real_,
TRUE ~ Collection_window_offspring)) %>%
filter(!is.na(Collection_window_offspring),
Female_mortality != 1 | is.na(Female_mortality))
# Create a function to build HTML searchable tables
my_data_table <- function(df){
datatable(
df, rownames=FALSE,
autoHideNavigation = TRUE,
extensions = c("Scroller", "Buttons"),
options = list(
dom = 'Bfrtip',
deferRender=TRUE,
scrollX=TRUE, scrollY=400,
scrollCollapse=TRUE,
buttons =
list('pageLength', 'colvis', 'csv', list(
extend = 'pdf',
pageSize = 'A4',
orientation = 'landscape',
filename = 'Productivity_data_clean_2')),
pageLength = 8500
)
)
}
my_data_table(Productivity_data_clean %>%
select(Mother_strain, Father_strain, Cross, subpopulation, Block, Week, Generation, Treatment, Vial_setup, Count_conditions, Female_offspring, Male_offspring, Collection_window_offspring, Pre_window_female_offspring, Pre_window_male_offspring, Pre_window_offspring, Total_female_offspring, Total_male_offspring, Total_offspring))
Column explanations
Mother strain: we measured extinction and productivity
of families derived from 36 different crosses. Crosses were only
conducted between strains that experienced the same
Evolution_treatment. This column shows the strain that the
founding female of the subpopulation was derived from.
Father strain: we measured extinction and productivity
of families derived from 36 different crosses. Crosses were only
conducted between strains that experienced the same evolution treatment.
This column shows the strain that the founding male of the subpopulation
was derived from.
Cross: an identifying ID for each of the 36 combinations
of Mother strain and Father strain.
subpopulation: lines of flies founded by single
inseminated females in generation zero. We generated the next generation
of each subpopulation by crossing a single full-sibling dyad.
Block: the experiment was run in 2 distinct blocks,
using flies separated by 9 generations. Block 1 ran for 29 generations,
whereas Block 2 ran for 20 generations.
Week: the week that the generation of each block was
conducted in. Note that blocks overlap between weeks 9-29.
Generation: the generation of the extinction assay,
ranging from 1 to 20.
Treatment: the strains had been exposed to one of three
evolutionary conditions for 20-29 generations: a female-limited response
to selection, a male-limited response and a control condition where an
evolutionary response occurred in both sexes.
Vial setup: specifies whether a vial was setup for a
given subpopulation in a given generation. Vials were not setup when the
subpopulation was extinct, or when it could no longer be included in the
experiment due to escape/ handling error.
Count_conditions: Normal indicates that the breeding
window lasted the standard three days. Extra day indicates the breeding
pair were left in the vial for an extra day. This occurred during
generation 15 of the second experimental block.
Female_offspring: the number of adult females produced
by the breeding pair that eclosed during the four-day collection
window
Male_offspring: as above, but for male offspring.
Collection_window_offspring: the sum of
Female_offspring and Male_offspring.
Pre_window_female_offspring: females progeny that
eclosed prior to the collection window opening.
Pre_window_male_offspring: as above, but for male
offspring.
Pre_window_offspring: the sum of
Pre_window_female_offspring and
Pre_window_male_offspring.
Total_female_offspring: the sum of
Female_offspring and
Pre_window_female_offspring.
Total_male_offspring: as above, but for male
offspring.
Total_offspring: the total number of offspring produced,
summed across the sexes and the pre-window and within-window
collections.
\(~\)
Do offspring emerge prior to day 10 throughout the experiment?
Productivity_data_clean %>%
ggplot(aes(x = Generation, y = Pre_window_offspring/Total_offspring)) +
geom_jitter(aes(size = Total_offspring), colour = met.brewer("OKeeffe2")[2], fill = met.brewer("OKeeffe2")[3], width = 0.2,
shape = 21, stroke =1, alpha = 0.25) +
scale_size(range = c(0, 5), name = "Total offspring") +
scale_y_continuous(limits = c(0,1), expand = expansion(mult = c(0.01, 0.01))) +
theme_tidybayes() +
coord_cartesian(xlim = c(1, 20)) +
labs(y = "Prop. offspring eclosing prior to day 10",
x = "Generation of inbreeding",
) +
theme(text = element_text(size = 14),
legend.position = "bottom")
Figure S3. Raw proportion of productivity that occurred prior to the eclosion window, for the first 20 generations of the experiment.
\(~\)
\(~\)
As mentioned above, we have strong expectations for how full-sibling
inbreeding should affect productivity. Decreasing genome-wide
heterozygosity and the expression of recessive deleterious alleles that
comes with this should follow a pattern of exponential decay. This
non-linear process can be modelled using the brms
non-linear syntax.
Here we model the productivity at generation \(t\) via the function \(a - e^{\lambda t}\) where \(a\) and \(\lambda\) are parameters that control the
initial productivity and the rate of its decay as generations of
inbreeding progress. We force \(a\) to
be positive using the exp() function, which expresses it on
the exponential scale. Large \(a\)
indicates higher starting productivity at generation 0 (we started
observation at generation 1, the productivity for which is given by
\(a - e^{\lambda}\)). Positive values
of \(\lambda\) indicate that
productivity declines with inbreeding, which is our strong expectation
and clearly the case following inspection of the raw data (see Figure
SX).
We model inbreeding depression while including zero values for all rows where a subpopulation was extinct. This helps combat underestimation of the decay parameter, caused by the non-random extinction of subpopulations, where high fitness subpopulations become over-represented as the generations of inbreeding progress (i.e. selection).
Fixed and random effects
In the brms non-linear syntax, different formulas can be
fit for each parameter. We fit the same fixed effect for \(a\) and \(\lambda\):
Evolution treatment. We fit subpopulation as a
random effect to account for repeated measurements taken on each
subpopulation for \(a\), but are unable
to do this for \(\lambda\), as
individual subpopulation productivity curves are highly heteroskedastic,
and are skewed by extinction events. To minimise pseudoreplication for
\(\lambda\), we fit Cross
as a random effect, as we do for the extinction data. We also fit
Week as a random effect for both parameters; this controls
for environmental variation between generations.
Priors
We fit moderately informative priors, with the aim to regularise our posterior estimates and improve model fitting by ruling out non-sensical values. Note that priors are particularly important when fitting non-linear models.
Note that \(\alpha\) related priors are expressed on the exponential scale.
\(\alpha\) (initial productivity) ~ Normal(\(\mu\) = 4.2, \(\sigma\) = 0.25)
\(\lambda\) (the rate of decay) ~ Gamma(shape = 0.8, rate = 4)
\(\sigma\) for \(\alpha\) ~ Gamma(shape = 2, rate = 1.5)
\(\sigma\) for \(\lambda\) ~ Gamma(shape = 0.25, rate = 2)
\(~\)
Prior predictive simulation for the productivity decline model. The plot shows the results from 100 prior draws.
n <- 1e4
# define the x-axis breaks
at <- 1:20
# how many prior draws would you like?
n_draws <- 200
# simulate
set.seed(16)
prior <-
tibble(index = 1:n,
a = rnorm(n, mean = 4.2, sd = 0.25),
lambda = rgamma(n, 0.8, 4)) %>%
slice_sample(n = n_draws) %>%
expand(nesting(index, a, lambda),
Generation = seq(from = 0, to = 20, length.out = 1e2)) %>%
mutate(productivity = exp(a - lambda * Generation))
prior %>%
ggplot(aes(x = Generation, y = productivity, group = index)) +
geom_line(size = 1/2, alpha = 1/4, colour = met.brewer("Hiroshige")[2]) +
theme_tidybayes() +
scale_y_continuous(limits = c(0,200), expand = expansion(mult = c(0, .1))) +
scale_x_continuous(limits = c(0,20), expand = c(0, 0)) +
labs(y = "subpopulation productivity", x = "Generation of inbreeding") +
theme(text = element_text(size = 14))
\(~\)
nonlinear_productivity <-
brm(bf(Collection_window_offspring ~ exp(a - (lambda * Generation)),
a ~ 1 + Treatment + (1|Week) + (1|subpopulation),
lambda ~ 1 + Treatment + (1|Week) + (1|Cross),
nl = TRUE),
data = Productivity_data_clean,
prior = c(prior(normal(4.2, 0.25), class = b, coef = Intercept, nlpar = "a"),
prior(normal(0, 0.25), class = b, coef = TreatmentFemale, nlpar = "a"),
prior(normal(0, 0.25), class = b, coef = TreatmentMale, nlpar = "a"),
prior(gamma(0.8, 4), class = b, coef = Intercept, nlpar = "lambda"),
prior(normal(0, 0.2), class = b, coef = TreatmentFemale, nlpar = "lambda"),
prior(normal(0, 0.2), class = b, coef = TreatmentMale, nlpar = "lambda"),
prior(gamma(2, 1.5), class = sd, nlpar = "a"),
prior(gamma(0.25, 2), class = sd, nlpar = "lambda")),
control = list(adapt_delta = 0.9, max_treedepth = 12), seed = 5,
iter = 20000, warmup = 5000, chains = 4, cores = 4,
file = "Fits/non_linear_productivity_model_no_fam")
\(~\)
\(~\)
new_data <-
Productivity_data_clean_2 %>%
select(Generation, Treatment) %>%
distinct() %>%
mutate(key = paste("V", 1:n(), sep = ""))
preds <-
as.data.frame(fitted(nonlinear_productivity, newdata = new_data, re_formula = NA, summary=F)) %>%
mutate(draw = 1:n()) %>%
gather(key, Collection_window_offspring, -draw) %>%
as_tibble() %>%
left_join(new_data, by = "key") %>%
select(-key) %>%
mutate(Treatment = case_when(
Treatment == "Female" ~ "Female-limited",
Treatment == "Male" ~ "Male-limited",
Treatment == "Control" ~ "Control"))
preds_summary <-
preds %>%
group_by(Treatment, Generation) %>%
median_qi(Collection_window_offspring, .width = 0.5)
nl_p1 <-
Productivity_data_clean_2 %>%
ggplot(aes(x = Generation, y = Collection_window_offspring)) +
geom_jitter(colour = met.brewer("OKeeffe2")[2], width = 0.25,
shape = 1.5, stroke =2, size = 1.2, alpha = 0.2) +
# geom_smooth(size = 1.5, colour = met.brewer("OKeeffe2")[7], alpha = 0.4) +
geom_ribbon(data = preds_summary,
aes(ymin = .lower, ymax = .upper, fill = Treatment), alpha = 1/2) +
geom_line(aes(group = Treatment, colour = Treatment), data = preds_summary,
size = 1, alpha = 1, linetype = 5) +
scale_linetype_manual(values = c(1, 2, 3)) +
scale_colour_manual(values = met.brewer("Hiroshige", 3), guide = "none") +
scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
scale_y_continuous(expand = expansion(mult = c(0.01, 0))) +
theme_tidybayes() +
coord_cartesian(xlim = c(1, 20), ylim = c(0, 150)) +
labs(y = "Subpopulation productivity", x = "Generation of inbreeding", fill = "Inheritance treatment") +
theme(text = element_text(size = 14),
legend.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.position = c(.99, .99),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6))
\(~\)
Estimating Generation one productivity and \(\lambda\)
Our model estimates \(a\) the initial productivity parameter. This is the productivity in generation 0, which in practice is non-sensical. We instead present productivity in generation one which = \(a - e^{\lambda*1}\).
Gen_1_productivity <-
preds %>%
filter(Generation == 1)
Initial_productivity_plot <-
Gen_1_productivity %>%
ggplot(aes(x = Collection_window_offspring, y = Treatment)) +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 1,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
labs(x= "Productivity in generation 1", y = "Inheritance treatment") +
theme_bw() +
theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
panel.grid.minor.x = element_blank(),
legend.position = "none", #transparent legend panel
text = element_text(size=14))
decay_estimates <-
as_draws_df(nonlinear_productivity, variable = "^b_lambda", regex = TRUE) %>%
mutate(`Female-limited` = b_lambda_Intercept + b_lambda_TreatmentFemale,
`Male-limited` = b_lambda_Intercept + b_lambda_TreatmentMale,
Control = b_lambda_Intercept) %>%
select(`Female-limited`, `Male-limited`, Control) %>%
mutate(posterior_sample = 1:n()) %>%
gather(Evolution_treatment, lambda, -posterior_sample)
decay_plot <-
decay_estimates %>%
ggplot(aes(x = lambda, y = Evolution_treatment)) +
stat_halfeye(aes(fill = Evolution_treatment), .width = c(0.66, 0.95), alpha = 1,
point_interval = "mean_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5, scale = 0.8) +
scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
labs(x=expression(paste("Decay rate parameter ", lambda,)), y = NULL) +
theme_bw() +
theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
panel.grid.minor.x = element_blank(),
legend.position = "none", #transparent legend panel
text = element_text(size=14))
Create plots d and e
Difference_gens_productivity <-
preds %>%
pivot_wider(values_from = "Collection_window_offspring", names_from = "Treatment") %>%
mutate(diff.mc = `Male-limited` - Control,
diff.mf = `Male-limited` - `Female-limited`,
diff.fc = `Female-limited` - Control) %>%
select(Generation, contains("diff")) %>%
pivot_longer(cols = contains("diff"), values_to = "productivity_diff", names_to = "Difference contrast")
# Make the three plots
Leaf_plot_mc_productivity <-
Difference_gens_productivity %>%
filter(`Difference contrast` == "diff.mc") %>%
ggplot(aes(productivity_diff, as.factor(Generation))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Purp") +
coord_cartesian(xlim = c(-5, 25)) +
scale_y_discrete(limits=rev) +
geom_vline(linetype = 2, xintercept = 0, size = 1) +
labs(x = "Male - Control\nProductivity diff.",
y = "Generation of inbreeding") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
Leaf_plot_mf_productivity <-
Difference_gens_productivity %>%
filter(`Difference contrast` == "diff.mf") %>%
ggplot(aes(productivity_diff, as.factor(Generation))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "TealGrn") +
coord_cartesian(xlim = c(-5, 25)) +
scale_y_discrete(limits=rev) +
geom_vline(linetype = 2, xintercept = 0, size = 1) +
labs(x = "Male - Female\nProductivity diff.",
y = NULL) +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
Leaf_plot_fc_productivity <-
Difference_gens_productivity %>%
filter(`Difference contrast` == "diff.fc") %>%
ggplot(aes(productivity_diff, as.factor(Generation))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Peach") +
coord_cartesian(xlim = c(-5, 25)) +
scale_y_discrete(limits=rev) +
geom_vline(linetype = 2, xintercept = 0, size = 1) +
labs(x = "Female - Control\nProductivity diff.",
y = NULL) +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
Combine the 6 panels
(nl_p1 / (Initial_productivity_plot + decay_plot)) / (Leaf_plot_mc_productivity + Leaf_plot_mf_productivity + Leaf_plot_fc_productivity) +
plot_annotation(tag_levels = 'a') +
plot_layout(heights = c(1.5, 1, 1))
Figure 2. Panel a shows productivity decline under increasing stress from inbreeding, split by inheritance treatment. Bands are 50% credible intervals to aid visual inspection and circular points show raw productivity counts for each sub-population in each generation. b shows the posterior distribution of productivity in the first generation of the experiment; white points are the median with associated 66 and 95% credible intervals. c shows the distribution of the productivity decay rate for each treatment. d-f show the difference contrasts for productivity in each generation. The innermost colour band approximates the median, while the middle and outermost bars show 66 and 95% credible intervals.